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Abstract 

We study various correlation functions (two and three point functions) in a large N matrix 
model of six commuting matrices with a numerical Monte Carlo algorithm. This is equivalent to a 
model of a gas of particles in six dimensions with a confining quadratic potential and logarithmic 
repulsions at finite temperature, where we are measuring the leading order non-gaussianities in the 
thermal fluctuations. This is a simplified model of the low energy dynamics of M = 4 SYM at 
strong coupling. We find strong evidence that the simplified matrix model matches with the dual 
gravitational description of three point functions in the AdS/CFT correspondence. 
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I. INTRODUCTION 

The AdS / CFT correspondence [l| in its simplest setting states that an ordinary quantum 
field theory in d — 4 dimensions, the maximally supersymmetric Yang Mills theory with 
gauge group U(N), is equivalent to type IIB string theory as a theory of gravity compactified 
on an AdS§ x S 5 geometry. As such, this equivalence provides a way to answer questions 
about the nature of quantum gravity by performing calculations in the dual quantum field 
theory. 

The radius R of the S 5 geometry in string units is related to the 't Hooft coupling of the 
field theory A = g 2 N via a scaling R 4 ~ A. The S 5 sphere is large in string theory units if 
R » 1. In this regime, it is expected that semiclassical calculations in supergravity are a 
reliable description of the physics. In other regimes, it is not known how to do calculations 
in the gravity side of the correspondence because one would need to understand how stringy 
corrections affect the dynamics of gravity. Thus, in order to compare the field theory and 
quantum gravity one should expand the field theory quantities at large values of A. This is, 
one should analyze the field theory at strong coupling. 

On general grounds studying field theory dynamics at strong coupling is a hard problem. 
In lattice gauge theory occasionally there are ways to address this problem Q], although 
since the M = 4 SYM theory does not confine (it is a non-trivial conformal field theory) a 
lattice definition of the field theory might not be very useful and it might be hard to extract 
information from such a formulation. 



A proposal for how 
a S 3 was proposed in 



;o do a strong coupling expansion of the Af = 4 SYM compactified on 
s|. The compactification on S 3 provides an infrared regulator of the 
dynamics. Morever, via the operator state correspondence, the spectrum of the Hamiltonian 
of the field theory on the S 3 is equivalent to the spectrum of dimensions of local operator 
insertions in Euclidean space. Thus, one can address the problem of computing operator 
dimensions and correlation functions by using Hamiltonian methods. 

In particular, in the proposal of [3j an expansion of the action of all fields of the M = 4 in 
spherical harmonics on the S 3 was truncated to the s-wave of the scalar field modes in the 
field theory. One also needs to include the s-wave of the Aq component of the gauge field to 
be able to impose the gauge constraint. The intuition that led to this truncation was that 
supersymmetric states in the free field limit only excite these modes and that supersymmetry 
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should account for cancelations of quantum corrections. Thus a naive truncation might be 
good enough to describe the dynamics, even at strong coupling. 

The effective Hamiltonian for this setup is a matrix quantum mechanics of six hermitian 
N x N matrices X, and their momenta 

where the coupling constant A has been calculated in 0]. It is proportional to the gauge 
coupling constant squared. This has to be supplemented by the gauge constraints. At large 
N, one expects that the eigenvalues of X are of order y/N. The first two terms of the 
Hamiltonian would be of order N 2 , but the last term, the commutator squared term, would 
be of order g 2 NN 2 = XN 2 , so for general configurations of matrices of this size the potential 
term is much larger than the kinetic term. Under these circumstances one expects that one 
should expand the system around configurations that minimize the large potential term. 
This is, one should expand around configurations of commuting matrices. The expansion 
that was proposed in [3| is exactly to solve the reduced (gauged) matrix quantum mechan- 
ics model of commuting matrices and to treat all other modes perturbatively around this 
(non-perturbative) background. Commuting matrices can be diagonalized simultaneously 
by a gauge transformation, so that only the eigenvalues are dynamical variables. For each 
diagonal component there are 6 coordinates, the eigenvalues of each matrix (we will label 
them as Xj). 

It was found that the wave function of the ground state for this reduced set of degrees of 
freedom was a Gaussian 

Vo = exp(-i^f 2 ) (2) 

i 

but that there is also a measure term that affects the calculation of averages in the quantum 
problem. This measure is 

= Y]_ Wi ~ Xj\ 2 ( 3 ) 

i<j 

a generalization of the Van der Monde determinant. 

Quantum averages of operators that are restricted to this set of diagonal variables are 
computed by the following integrals 

= J> 2 exp(-E^ 2 )Q(x) = / exp(-£ t x 2 - log |fl ~ S^)Q{x) 
jVexp(-Ei^) jV ex P(~ Ei ^? ) 



This is equivalent to studying thermal correlations for a Boltzman gas of N particles in six 
dimensions in the presence of a one body potential (x 2 ), and with a two body logarithmic 
repulsive interaction (\og\x — y\ 2 ))- This is a somewhat unusual problem in statistical 
mechanics, but it can be approached numerically. This approach was initiated in |5|, where 
very few quantities were computed, confirming some of the theoretical ansatz [3|, |4| that 
solved the saddle point approximation for the thermodynamic limit of the distribution (see 
also [6| for some related problems in the thermal field theory case). In particular, we expect 
that in the thermodynamic limit the theory will be described by some density of particles 
p(x) and thermal fluctuations of the density. 

In this paper, the numerical study of this ensemble is continued. In particular we compute 
the density fluctuations of various modes of the ensemble for various values of N and we 
study the approach of the fluctuations to the thermodynamic limit. We also compute the 
leading order non-gaussianities of the density fluctuations. These are related to three point 
(extremal) correlation functions in the conformal field theory. These have been studied in the 
general case |7| where a non-renormalization theorem for BPS operators was conjectured. 
However, it was found that in the extremal case the computation requires solving some 
subtleties [8[ and once this is done, the results of |7| for supergravity can be used. A non- 
renormalization theorem for these extremal correlators was proved in [9( using the techniques 
of harmonic superspace. 

The point of view we will take in this paper is that we want to verify if the approximations 
that were taken to study the dynamics of Af = 4 at strong coupling in [3] are a complete 
description of the low energy dynamics that leads to supergravity: we do not need the other 
modes in the field theory to reproduce the results expected by the non-renormalization 
theorem at strong coupling, and we can compare with the results found using supergravity 
in 0,3- 

The paper is organized as follows. In section two we explain how the two and three 
point functions are calculated in the free field theory limit of M = 4 SYM. In particular, we 
show how the OPE coefficients are computed by a free matrix model. We also discuss in this 
section how this calculations was matched to supergravity, suggesting a non-renormalization 
theorem for three point functions. Next, in section [TTT1 we give a description of how we set 
up the calculations of the three point functions for strong coupling within the wave function 
approach. In section IIVI we give a detailed presentation of our numerical results for two 



4 



and three point functions at strong coupling. Then we conclude. We have also included an 
appendix where we describe how the statistical error bars were calculated and some details 
of how the Monte-Carlo code was calibrated 1 . 

II. FREE MATRIX MODEL RESULTS 

Conformal field theories are usually determined by the two and three point functions of 
primary operators (these are given by local insertions of composite operator Oi(x)). The 
two point functions (of scalar operators) are given by 

(OiixWv)) ~ i (5) 

where Ajj are the dimensions of the corresponding operators, while CV,- is a set of numbers 
(symmetric in and the three point functions of scalar operators are given by 

(OiWOiivW)) ~ | x _ y |A <+ A J -A fc | :E _^A fc -A j | y _ g |A t+ A J -A < ( 6 ) 

where the are structure constants of the theory. 

The matrix CV,- is called the Zamolodchikov metric, and it is a standard procedure to 
diagonalize the metric so that = 5 it j. In this sense the exact normalization of the two 
point functions is somewhat unphysical, unless there are degeneracies in the list of operators 
with given quantum numbers. 

The coefficients are also the coefficients of the Operator product expansion of Oi 
and Oj into Ok (or any permutation of the three). 

O^OM-Y. u.a£%-a» 0*(*) (7) 
k 1 1 

This is usually an asymptotic series, organized starting from the most singular terms 9those 
with smallest value of A&. 

In the special case where for the leading order term A& = A, + Aj, the operator product 
expansion is non-singular, and one can take the limit x — > in the above formula. This 
makes it possible to have a composite operator (OiOj)(x). 

1 The computer code used to generate the data with intstructions for compilation is available on request 
from D. B. 
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In the special case of Af = 4 SYM, there is a special class of operators of protected 
dimension. This is, the dimension of the operator when calculated in the free field theory 
limit is identical to the dimension of the operator in the interacting theory. 

The simplest such operators are given by symmetric traceless tensors of 5*0(6), 

0, t , - Tr(X ll ...X / J(0) (8) 

where <Y'0;.j ,-, = and = C) M ...., ,. . This is a composite field of k free 

fields, and it is of dimension k. These were first described in the AdS/CFT setup in flo| . 
where they were matched with supergravity fluctuations Each trace in interpreted as 
a single graviton state. 

Each of these single trace operators can be related to a simpler one Tr(Z k ), where we 
choose a particular complex combination of the fundamental scalar fields Z = X\ + 1X2. 
This can be thought of as the highest weight state of the associated SO (6) representation. 

The (free field) two point functions of these operators depend on how we choose to 
normalize the fields in the field theory action. 

In free field theory, to leading order in planar diagrams, one finds that 

( Tr(Z k )(x) Tr(Z k )) ~ _(1 + 0(l/iV 2 )) = ^ (9) 
More precisely, one can compute C^k by performing integrals in a Gaussian matrix model 



for a complex matrix z 12 1 



where the average in the right hand side is understood as a matrix model integrals (overlap). 

One also finds that multi-traces of Z are protected operators, and they are different than 
single trace operators. Thus one has the problem that there is more than one operator of a 
given dimension for a given set of quantum numbers. In general the operators mix and this 



problem needs to be resolved. A complete solution in the free field limit was given in [121 ]. 
where a complete orthogonal basis was found. The general problem of multitrace mixing 
can also be understood in terms of the complex matrix model, by calculating averages of 
the form 

j{dzdz) N2 exp(- Tr(zz))Uj 1i:(z*0 II,-' Tr (^ 
j{dzdz) N2 exp(- Tr{zz)) 

(Q Tr(^)Il Tr ^')> ( 12 ) 
j 3' 



If ki + • • • + k s is held fixed, the numbers C provide a non-diagonal metric on the space of 
half BPS states. 

Of particular interest to us is the correlation function C n m „ n;m , and to leading order this 
is given by 

Cn,m-n;m = <m - ^ruN™" 1 (1 + 0(1/N 2 )) (13) 

in the same normalization for the fields where the two point functions are calculated by 

C k , k = kN k (l + 0(l/N 2 )) (14) 

C n ,m-n.-m can be interpreted both as a two point function for describing mixing between 
the operators Tr(Z n ) Tr(Z m ~ n and Tr(Z m ), and as a three point function for extremal 
correlators. 

To study the physically relevant information encoded in these two and three point func- 
tions it is natural to calculated them in a case where the two point functions are normalized 
relative to the Zamolodchikoiv metric for two point functions |7j. Thus, we would compute 
the normalized correlator 



N n>m -n;m = J= = ~ ^ ^ (1 + 0(l/iV 2 )) (15) 

The non-renormalization conjecture of [7J is that iV n>m;m _ n does not change in value in 
extrapolating from free field theory to strong coupling. This is the set of numbers we will 
compute at strong coupling with the proposal jjj. 

The other useful calculation to do is the following 

0„,...,, : ,,,„.,, ~ C s Kk s\{\ + 0(1/N 2 )) (16) 

where we have take s copies of the same trace. It is easy to show by planar counting that 
the right hand side in the above equation has a factor of s\. This is because the leading 
contractions are from disconnected diagrams between the different groups of traces. The 
factor of s! just counts how the different permutations of pairings of the traces. 

The equation ( TT6l) states that to leading order in 1/N, Tr(z k ) can be interpreted as a 
raising operator in a harmonic oscillator Tr(z k ) ~ a k . This is because 

ll(«Ir|0>|| 2 ~ S ! (17) 
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thus if we think of equation (fTOj) as a statistical ensemble, we should find that the distribution 
of Tr(z k ) is Gaussian. This is familiar from quantum mechanics for a harmonic oscillator, 
where the ground state wave function is 

i/>o(x) ~ exp(-x 2 /2a) (18) 

and we find that the distribution of x in the ground state is Gaussian. This also applies 
to the combination a* ~ x + ip, where we now get a Gaussian distribution in the complex 
plane (phase space for a single variable). 

One point functions of powers of a complex matrix vanish 

( Tr(^)> = (19) 

for n > 0. This is because the gaussian measure is invariant under phase shifts z — > exp(i0)z. 
This property is inherited by all correlation functions, so correlation functions that are not 
invariant under that shift will automatically vanish. 



III. STRONG COUPLING SETUP 



We want to give a prescription to calculate the three point functions for the M = 4 SYM 

n 

at strong coupling within the proposal of [3|]. The idea is to exploit the relation between the 
Gaussian matrix model and the OPE coefficient. 

We should notice that the formulae for computing overlaps in equations (IIOIIIII) is very 
similar to the formulae for computing averages in the thermal gas in equation (j3j). 

Indeed, the main setup of states that half BPS wavefunctions associated to operators 
Tr(Z n ) are computed at strong coupling by taking the ground state wave function in di- 
agonal variables and multiplying it by Tr(Z n ) where Z = X\ + 1X2 is a complex diagonal 
matrix (or other such product of traces). Also, the gaussian matrix measure is associated to 



13|, 



the ground state wave function for the s-wave of the scalar fields in the free field limit 
when the field theory is compactified on a 5 3 . 

This same interpretation for the measure H] is available. Thus it is natural to compare 
the vacuum expectation values of equations (BJ and ( TTUl fTTTl . 

Under this comparisson the same type of formula as above holds, where we substitute 

Tr(Z") = 5>J + i^r~ f z n p (20) 



This is, the traces are calculated in the matrix model of commuting matrices by summing 
over the eigenvalues. They can also be evaluated by calculating integrals over the density 
of particles in the Boltzman gas in the thermodynamics limit. Thus the integrals represent 
moments of the density distribution. 

The density distribution of particles is given in the saddle point (thermodynamic limit) 
by flfl 

p (x) = (p(x))~5(\x\-r ) (21) 

where ro = \J N/2, and the density is normalized so that J p = N. 

We can think of the true density as given by p = po + 8p, where 5p are density fluctuations 
and (8p) = 0. 

It is also convenient to normalize the size of the sphere so that the value of vq is scaled 
to one. This is done by rescaling x = x/r . In this way we find that 

Tr(z n ) = J p(x)z n ~ J p(x)Y n (6) (22) 

This is, the normalized traces identify the integral of the density times a particular spherical 
harmonic of the sphere S 5 . As such, these traces compute angular fluctuations of the particle 
density on some specific spherical harmonics. We will use the notation Sp(n) to indicate 
schematically the density fluctuations integrated over the corresponding spherical harmonic, 
so that 

J p(x)Y n (6) ~ 8p(n) (23) 
The idea is to study the two point density fluctuations 



p{x)Y n {9) I p(x)Y:(9)) ~ {5p\n)) (24) 
and the three point density fluctuations 

( / p(x)Y n (6) [ p{x)Y m _ n {9) / p(x)Y; n (6)) ~ (5p(n)5p(m - n)8p{m)) (20) 



The first set of quadratic fluctuations can be considered the Gaussian fluctuations: given 
the numbers (5 2 p(n)) there is a unique Gaussian with central value zero whose fluctuations 
are characterized by those numbers. 

The three point functions would vanish if the spectrum is fully Gaussian between the 
traces. Thus a non-vanishing of the three-point functions reflects a failure of the density 
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distributions to be Gaussian. Since we expect the normalized three point functions to be of 
order 1/N (this follows from comparing to equation ffl5l) . and from standard 1/N counting 
arguments [ijj]), we expect the non-gaussianities to be small in the large N limit. The 
leading order non-gaussianities are these three point functions. 

The method for calculating the three point functions is to generate a set of configurations 
distributed according to the measure H] via a Monte Carlo algorithm. The average value of 
the two and three point functions can be calculated by averaging over a large sample of 
statistically independent typical configurations generated by the Monte-Carlo method. One 
can also sample the distribution of the 5p(n) directly and check that it is Gaussian. The 
algorithm used was presented in [5]. Here we use the same algorithm for various values of 
N. This will let us study the approach to the thermodynamic limit, to estimate the 1/N 2 
corrections. It is also useful to point out that a similar ensemble was studied in [15], where 

n 

a similar measure term was obtained. However, in that paper [15J there is no gaussian 
confining potential for the eigenvalues. This gaussian factor simplifies the analysis of the 
large N limit considerably. 



IV. NUMERICAL RESULTS 



We wish to study typical configurations of the gas of particles described by the equation 
(J1J), and to evaluate the corresponding distributions for the possible values of various 0(x). 
We do this by studying individual configurations of the gas, with a measure proportional to 

dfi ~ exp(- + ^log \$i ~ Xj\ 2 ) (26) 

p i<3 

The Monte-Carlo algorithm to navigate these configurations was described in Jsj] . We move 
from one configuration to the next one by applying a random move on particle i, where each 
of the coordinates of Xi is varied by a random number between [—0.55,0.55]. We accept a 
new configuration according to the Metropolis-Hastings criterion. The typical values of 5 
that we use are between 4 and 7. We do this by cycling over all the particles in the gas and 
updating them one at a time. Roughly 40% of the particles are updated in each cycle. 

We found that storing data from the configurations every 30 cycles was optimal. By 
this time the autocorrelations of the ensemble have died down considerably and for all 
practical purposes the configurations can be considered to be statistically independent. This 

10 



is discussed in more detail in the appendix. 

To obtain more information from each configuration, we chose various different orienta- 
tions of the complex coordinate z, where 

z = x a + ix 13 for a < f3 (27) 

This choice is a symmetry of the system and it allows us to gain efficiency in the computation 
of the correlation functions. We wrote down only the multipoles themselves and not the full 
configuration of particles. 

For two point functions, we did the calculations with iV = 400, IK = 
10 3 , 2K, 4K, 10K, 25K particles. For three point functions, which require more statistically 
independent configurations to measure numbers statistically different than zero, we were 
able to get as high as 5K particles. We studied the first ten multipole coefficients Sp(i), for 
% = 1,...10. 

A. Two point functions 

The first measurements we made on the sample where those of the two point functions. 
These are expected to have a Gaussian distribution for Sp(n) ~ J 5pY n (6). The reason for 
this is that the large N counting responsible for equation [16] is robust independent of the the 
value of the 't Hooft coupling. We expect to find a Gaussian distribution for each 5p{n) in 
both the real and imaginary part. Because of the invariance of the ensemble under rotations 
of z by phases, the real and imaginary part distribution will be the same. Also, the mean 
of the distribution vanishes. 

To test for Gaussianity, we plot a histogram of the distribution in the interval [—8a, 8a] 
for the combined real part and imaginary part, divided into fifty intervals. This is depicted 
in figure [1] for our largest statistical sample. 

One can put all ten multipole histograms on top of eah other, and they are statistically 
indistinguishable from one another. One also finds that the data is very well described 
by a Gaussian distribution. Indeed, it must be the case that the distribution for Sp(l) is 
Gaussian. This observation follows from writing the partition sum in terms of a center of 
mass variable xqm = N' 1 J2 P %pi an d relative coordinates. The two particle potentials are 
independent of xqm- The quadratic function J2 p ^p can be written as a sum of squares of 
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FIG. 1: Sample histogram distributions for 5p(n), N=5K. n = 1, 2 

relative distances plus a quadratic part. Thus the distribution for Sp(l) is Gaussian. This 
is also familiar from the full field theory for U(N): the diagonal U(l) degrees of freedom 
decouple completely. 

If we normalize the radius of the distribution to one, we find that the fluctuations in the 
center of mass coordinate have a theoretical Standard deviation identically equal to one. We 
use this as a consistency check on the performance of the code. The fact that the histograms 
for all multipole moments look identical means we have a high degree of confidence that all 
the distributions are gaussian. 

We can also look at the two point function absolute normalization, and how it depends 
on N. The results we find are depicted in figure [3 Statistical error bars are small and are 
not shown. It is computationally cheap to generate large data samples to evaluate these 
precisely. Here we are interested in the qualitative approach to the large N limit, not on 
measuring deviations from the large N limit precisely (we expect that these should be a 
series in 1/N 2 and are interesting on their own right, as they measure some type of quantum 
corrections in the dual theory, but their study is beyond the scope of the present paper). 

The results for the absolute normalization of the two point functions are surprising in 
that the approach to the large N limit for the individual multipole coefficients is slow. One 
can be off by 20% in the normalization of <5(10) 2 for iV = 10 4 . This would suggest that 1/N 2 
corrections are rather large, even though N is very large. 

Indeed, the pattern seen in the figure [2] suggests that the two point functions have a local 
maximum that depends on N. This local maximum moves to higher multipole moments as 
we increase N. Moreover, the maximum seems to occur at about half the large N limit value 
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FIG. 2: Two point functions for N=25K, N=10K, N=4K, N=2K, N=1K, N=400 respectively, from 
top to bottom. Labeled by angular momentum label n 

for the corresponding multipole moment (in our simulation the large N limit is declared to 
be N = 25K), so it can serve as a proxy for the scaling of the 1/N corrections. We find that 
the value of J where the local maximum of S( J) 2 occurs seems to scale like iV 1 / 4 (the best fit 
to the data gives us J ~ jV 1 ' 3 - 9 , but the value of 3.9 has moderate systematic errors because 
J can only take discrete values and we are not finding sufficiently large values of J where 
these errors could become negligible). This is, corrections are suppressed by powers of J A /N 
(or equivalently J/iV 1 / 4 . The value of iV 1//4 found here matches nicely with the expectation 
from gravity. After all, the radius of the sphere in AdS§ x S 5 in the dual gravitational theory 
scales like iV 1//4 in Planck units. 

What this means is that corrections in J (which is interpreted as momentum on the 
sphere) are suppressed by the Planck scale. This is in essence another measurement of the 
Planck scale than the one found in Q], which was related to the back-reaction of the geometry 
to the presence of extended objects. This type of coincidences make it very plausible to 
believe that the matrix model of commuting matrices found in Q| is actually capturing all 
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the relevant physics of the strong coupling system that is described by a dual geometry. 
This is in contrast to the weak coupling (free field theory regime), where corrections scale 
like powers of J 2 /N 16 1. 

We should also contrast these results with the leading order two point function, as ex- 
pressed in equation 0, where it is seen that the two point function scales as ~ k, once 
the factors of N are removed by the normalization of the Z. Here we find that at strong 
coupling the value of this quantity scales differently with k. Indeed, the figure [2] suggests 
that the scaling is power-law. A fit to the data suggests that 



a 



k,k 



k 



1.75 



(28) 



more precisely, we make a fit for the power-law for various values of N, finding the following 
values, depicted in tabled] 



N 


400 


IK 


2K 


4K 


19K 


25K 


Exponent 


1.3 


1.48 


1.57 


1.61 


1.70 


1.74 



TABLE I: Fit to power law behavior Ckk ~ k a for different values of N 



The exponent 1.75 is guessed because it is a simple rational number close to the observed 
value 1.74 1.75 = 2 — 1/4. The factor of 1/4 in the denominator is also suggestive from 
the fact that the sphere scales like N 1 ^ in Planck units in the gravity dual. It would be 
interesting if one can prove mathematically that the matrix model of commuting matrices 
predicts this value (this is, one would partially solve the matrix model). 

B. Three point functions 

Now we turn to the problem of evaluating three point functions numerically. What we 
need to do is compute the average values of 

( Tr(i n ) Tr(z m " n ) Tr(l m )) ~ 0(1/ N) (29) 

and we expect the results to scale like 1/N. For an individual configuration, the typical 
value of Tr(z n ) ~ 5p{n) is found to be of order 1. This means that we need to measure 
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the center value of the corresponding distribution by a factor oll/N better than the typical 
width of the distribution. 

If we have a statistically independent sample of L draws of the distribution, we find that 
the central value of the the distribution is measured to accuracy 1 / y/L relative to the width 
(this is the central limit theorem). We need this precision to be of order 1/N to distinguish 
the central value of the distribution from zero. This is, in order to measure a non-gaussian 
behavior in the distribution we need large statistical samples. In our case L scales like N 2 . 
We should also notice that in the Monte-Carlo code, each loop for updating the particles 
require 0(N 2 ) computations. The CPU time required to get to the same precision in the 
measurement of non-gaussianities for higher values of N scales like iV 4 . 

The possibility of a non-renormalization theorem as described in {^J lets us compare with 
the values of these non-gaussianities in free field theory, as expressed by [151 This is, we 
can measure in the Monte- Carlo and compare absolutely to the result found in gravity by 
rescaling the three point function we measure by the expected value. Thus, we can compare 
universally for various values of N if we are near the expected value or not. Seeing as the 
two point functions behave differently at strong coupling than at weak coupling in figure 
|2l a match of the three point functions would be a strong test that the matrix model of 
commuting matrices is correct. Our results are presented in figure [31 

The most important number for our considerations is L, the number of statistically inde- 
pendent configurations. We have that for N = 5K we have a total of 242151 measurements 
with very small correlations. This is about a tenth of (5000) 2 , which is the suggested num- 
ber of configurations that we would need. This sample took over a month of CPU time 
to generate. What makes the calculation possible is that the coefficient we are trying to 
measure is of order \Jn{m — n)m larger than just 1/N. This number can be as large as 15, 
and is typically of order 7, so for these values we can cut the number of samples by 1/7 2 to 
get reasonable results. Thus, with the number of samples we have it is possible to test the 
commuting matrix model. 

The value that we are comparing to for each triplet of mutipoles is normalized to one. We 
plot the values of iV n m _ n;m that are measured in the statistical sample with a la error bar. 
We notice that the error bars become smaller as we increase m. This is expected because 
the corresponding three point functions are larger relative to the width of the distribution. 

On the other hand, we expect that as we go to higher multipole moments the "quantum 
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FIG. 3: Normalized three point functions, a) N= IK, b) N= 2K, c) N= 5K, with la error bars. 
The three integer labels in the x axis are the values of m, n, m — n multipoles that we are comparing. 

corrections" become larger and larger and we are further away from the large N limit value. 
We should also have in mind that the center of mass mode decouples completely (this is 
referred to as the singleton sector in the supegravity theory), so one should ignore the 
correlations with n — 1 for the most part. They are artificial in that we have defined Tr(z n ) 
without substracting the center of mass mode (a discussion of these issues can also be found 



in 



17|]). Indeed, this is what one would do if one were studying the SU(N) theory rather 



than U(N). In general, this is a 1/N correction to the traces, and these small corrections 
will be of order 1/iV 2 in the definition of the other three point functions. Thus we tabulate 
them in the unsubstracted form, which is computationally simpler, without making any 
significant errors. We left the data in the graph with n = 1 to show all of our numerical 
results. 

We see from the figure [3] that as iV increases, the match to the theoretical supergravity 
value improves. We do not have a clear understanding of the systematic deviations from 
being at finite N rather than in the thermodynamic limit (we will call this our systematic 
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error). In practice, we would need larger values of N and much larger statistical samples to 
be able to study this question numerically. 

However, we can estimate the size of the expected variations by using the deviations of 
the two point functions as given in figure [2] as a proxy for the systematic errors. This is, we 
take the non-normalized three point functions and divide them by the correct value of the 
two point functions in the thermodynamic limit (four our purposes this is N = 2bK. The 
deviations in the central value are depicted in figure HI The large value corresponds to the 
measured two point functions at the given value of iV while the smalll value corresponds 
to the two point function normalization for large N. A better theoretical estimate would 
have us fit the data to an expansion in 1/N 2 and extrapolate to iV — > oo. For this to 
work, we need to ensure that we are in a regime where only the 1/N 2 correction matters. 
Unfortunately, the two point functions seem to indicate that higher corrections in 1/N 2 are 
relevant. Also, our data quality is too poor to do this reliably and we consider the estimate 
described above more meaningful. 

----- : ■ ■ J ■ ' 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

o) : 



in, ,iii||l|l 

b) - 



D 2 

E , 

O 

z 

-1 

-2 

r-" n n n <t in 4 w id in •* n id* in ^" o" n" a) m" w, e n d_ w 

N ^ 4 in in ID IO Id" N N N ID ID CO 00 O) a" O) O) o o o o o 





A 




3 




2 


o 


'o 


1 









o 




o 


-1 


o 


-2 



o 

Id 
o 

O 

o 



. ■ . 11 ■ J ■ ■ ,,11 

c) : 



FIG. 4: A naive estimate of the systematic corrections for the same data of figure El 

From the figure H] we see that a naive analysis of the systematic corrections suggests that 
the values will move closer to the gravity value, and that these are becoming smaller as N 
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grows. 

Just for comparison, we include here a graph for the non-normalized three point functions 
in figure El for the particular case of N = 5K. As the graph shows, the pattern of three 
point functions has a lot of variation in the numbers that we measure (the range is over 
two orders of magnitude). The fact that in figure [3] the data becomes very flat and of order 
one suggests very strongly that modeling the data by the gravity prediction gives a very 
reasonable fit to the numerical data. 



1.5 



1.0 



0.5 



0.04 



I 5 i J 



FIG. 5: Non-normalized data for three point functions for N=5000 



Our final numerical results are that the matrix model calculation is consistent with the 
gravity computation. Since the systematic errors are still large one can not claim at this 
point that there is an exact match, but there is strong evidence for agreement. Considering 
that there were large corrections in the two point functions, we believe this is a strong test 
of the commuting matrix model proposal. 
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V. CONCLUSION 



In this paper we have studied the extremal three point functions in the commuting matrix 
model proposal for the strong coupling expansion of J\f = 4 SYM at strong coupling [3!]. We 
did this by calculating a statistical average of the ground state wave function of the model, 
following a Monte Carlo algorithm that was used in the previous work Q]. Considering 
that the values of the two point functions that we measured were very different than what 
was expected from the weak coupling calculations (they have different scaling with the 
multipole number), it seems rather surprising that the normalized three point functions 
match. Although this is expected from non-renormalization theorems 9J, it is hard to 
imagine that with all the approximations that are made to derive the model in [3J that the 
normalized three point functions are essentially uncorrected. 

We have also seen numerically that the 1/N corrections have different parametric de- 
pendence on the quantum numbers of a state than at weak 't Hooft coupling. This fact 
motivates us to try to understand better the proposal given in [3J. Recent work in studying 
more general examples of the AdS/CFT correspondence [3] has shown that this type of 
proposal can be extended to many other setups. It has already been shown that a big part of 
the supergravity spectrum including many non-BPS multiplets is matched. That is exactly 
the collection of states that we are able to simulate in theA/" = 4 SYM setup. 

Thus, in principle, if one is able to numerically evaluate the distributions of particles 
that generate other geometries, one would be in a position to study three point functions of 
primary operators on strongly coupled field theories numerically. Here one can do a strong 
test of the AdS/CFT correspondence by matching the three point functions of field theory to 
those of gravity. Seeing as in many of these field theories there is no perturbative limit where 
some of these calculations could be performed, it would be the first instance of calculating 
the field theory correlators (and therefore the OPE expansion coefficients) directly at strong 
coupling. To do this one needs some detailed information about the metric of Sasaki- Einstein 
manifolds, and it might be necessary to compute these numerically along the lines of 19] . 

It would also be interesting to obtain better statistics and a systematic determination of 
how the 1/N corrections behave for both two and three point functions in M = 4 SYM. 
These encode non-trivial higher genus corrections in the dual string theory. 
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APPENDIX A: STATISTICAL ERROR BARS AND CORRELATIONS 



In this appendix we give for completeness a brief description of the computation of the 
statistical error bars in our results, as well as how we optimized the extraction of data from 
the Monte Carlo generator. This is important because some of the measurements we are 
making (the calculation of three point functions) require calculating an average of a quantity 
to precisions much higher than it's typical fluctuation. Indeed, we are trying to measure 
quantities of order 1/JV to fluctuations of order one. 

The main observation is that when we compute an average in a Monte Carlo setup, we 
are computing the average of a time series, where the time in question is simulation time. 

Thus, we are computing 

t=l 

where L is the number of recorded loops. 

We want to estimate the statistical error bar for h. If the h(t) are statistically indepen- 
dent, then the variances satisfy 

a\(h)) ~ a 2 (h)/L (A2) 

However, in practice, the h(t) are correlated. The correlations can be computed by the 
covariance of the time series 

(6h(t + i)6h(i)) 

c(t) — W) — {) 

U sing c(t), one can compute the error in equation (1A1I) by using the standard formula 
(see [2Qj, eq. 5.13, for example) 

a 2 ((h))~L-^(h)(l + 2j2c(t)) (A4) 

t=i 

where n is a suitable cutoff. 
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Typically the c(t) decay exponentially, c(t) ~ exp(— At), where A depends on the observ- 
able. To some extent A -1 measures the relaxation time of the corresponding variable, so it 
is convenient to plot log(|c(i)|) with respect to simulation time for various variables. 




FIG. 6: Logarithm of Correlations of two point functions, for N=25K, steps per recording=5 

The figure [6] shows the logarithm of the autocorrelations of the different multipoles for 
the simulation with N = 25 x 10 3 particles, and recording the data after five simulation 
steps. The different curves represent the different multipoles. It is clearly seen that at the 
beginning there is a typical linear behavior, with different values of A for each multipole. 
There is usually a value where c(t) becomes oscillatory. This is a reflection of the limited 
statistical sample and it is a good place to cutoff the value of n used in the sum above. 

For the data shown in figure [3] the step used between recordings was 30 Monte Carlo 
iterations, where the two point autocorrelations decay very quickly (six times faster than in 
figure E]). We tuned the steps for that simulation so that roughly log(\c(l)\) < —1 for all 
two point function observables when we are doing three point function measurements. The 
statistical autocorrelation was used to write the error bars in figure [31 This resulted in a 
correction of less than about 3% for all values, if one assumed that the all the data taken 
was statistically independent. 

Although individual values of A depend on how one generates configurations, the ratios 
Xi/Xj might show a more stable pattern when one changes the parameters of a simulation. In 
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our case we see that the highest multipoles have longer autocorrelations. If one considers the 
A -1 numerically, they seem to show a pattern that scales simply with the mutipole number. 
This might be interesting to explore for the future. 
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